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Abstract 

We show that the ratio of different structure functions have a maximum which depends 
on xb and A. We argue that these maxima are proportional to the saturation scale. The 
analysis of leading and higher twist contributions for different observables is given with 
the aim of determining the kinematic region where high parton density effects could be 
seen experimentally. 
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1 Introduction 



It was noted many years ago that at fixed Q 2 but low xb (high energies) we are dealing with 
a parton system in DIS in a kinematical region where the coupling constant is small while 
the interaction between constituents is large 0. In such a parton system we have competition 
between emission of partons (gluons) and their annihilation, which leads to an equilibrium state 
with a definite value of the gluon density @||. At any value of small xb there is a value of 
Q 2 = Q 2 (xb) at which the parton density reaches a sufficiently high value that the interactions 
(annihilations) between partons become so strong that they diminish the number of partons. 
The new scale in DIS at which this phenomenon occurs we shall call the saturation scale. This 
scale depends on the energy of the colliding system as well as on the atomic numbers of colliding 
nuclei, and can be evaluated from the condition that the packing factor (f pa ck) of partons in a 
parton cascade is equal to unity [0, @, [3| : 

3ir 2 a s x b Ga(xb,Q 2 s (x b ;A)) 

J pack 0/02/' A\ D2 ' V ' 



2Q 2 s (x B ;A) ttR 



and r saturation = V^^b); 

The high parton density effects appear as nonlinear corrections to the DGLAP evolution 
equation. A few theoretical approaches have been developed which take into account non- 
linear corrections to the DGLAP evolution. One of them is based on the pQCD approach 
which has been used at the limit of the region where it is valid [[I], |2|, |J . Finally, Yu. Kovchegov 
and la. Balitsky |6] have derived an evolution equation for the dipole wave function with 
GLR-type corrections [l] which describes the parton system in the full kinematic region. The 
initial condition for this equation is a Glauber-like distribution which must be fixed at not too 
small values of xb- The saturation scale which appears in this equation has been studied in 
Refs.0. The same scale appears in a different approach based on the effective Lagrangian for 
semiclassical gluonic fields, emission and the interactions which are responsible for high gluon 
density system || ||. Despite the quite different techniques that have been incorporated in 
these approaches, they describe the same physics just approaching it from a different value of 
the parton density: pQCD from the low density side while the effective Lagrangian approach 
is from the high density side. 

Approaching the saturation scale the total cross section saturates at a geometrical size of 
the target. Using the language of the operator product expansion we can say that the reason 
for this breakdown is that it is no longer a valid approximation to neglect higher twist operators 
in comparison with the twist-2 operator. Thus, to estimate the kinematical region for which 
shadowing corrections , i.e. non-linear effects at large energies become important, one has to 
solve the non-linear evolution equation. 

This problem has not yet been solved, to do so we suggest the following approach. Namely, 
we propose estimating the size of the different twist contributions to the observed quantities and 
find the kinematical region where they become equally important. This will aid in determining 
the region where the linear evolution breaks down. Such an approach has been developed in 
Ref. for the DIS on a nucleon in the framework of Golec-Biernat and Wusthoff model |10| 



It should be stressed that in spite of the fact that at present there does not appear to be 
any direct evidence of saturation in the experimental data at HERA, there are some hints of 
high parton density effects in the same data [p]| . One of the goals of this paper is to provide 
experimentalists with new observables or/and with new ways of analyzing the data which will 
enable one to explore high density QCD effects in data from nuclear targets. 

The paper is organized as follows: We start by arguing that the eikonal approach is a 
reasonable model, since it is a first iteration of the exact evolution equation. We then derive, 
making an approximation, the analytical expression for the contribution of any twist to the 
longitudinal and transverse structure functions, as well as to the longitudinal and transverse 
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diffractive structure functions. We show that to estimate the saturation scale it is helpful to 
consider the ratios of the various structure functions. (Sec. 2). In Sec. 3 we discuss numerical 
estimations with formulae of Sec. 2. Finally, in Sec. 4 we conclude by proposing how one can 
measure different twists at eRHIC. 



2 The Model 

2.1 Eikonal Approach 

Complete description of the eikonal (or Mueller-Glauber) approach to shadowing corrections 
was given in Refs. [O, O, [TJ], [TR ITS]. We recall the main points. 



1. We work in the target rest frame where the whole evolution occurs in the virtual photon 
wave function. This wave function consists of colour dipoles which emerged as a result 
of numerous splittings of gluons (i.e. quarks and anti-quarks in N c ^> 1). In the leading 
log(l/xB) approximation a long time elapses after the photon wave function is formed, 
until every dipole scatters off the nucleus. The factorization of the two processes - 
radiation of the dipoles and their interaction with the target is the basis of the dipole 
picture of DIS[|l^, |i~3f . In this picture the scattering amplitude is diagonal with respect 
to the size of the dipole r t and the fraction of initial energy z carried by quark [Tj| . Let 



P 7 * (z, r t ; Q 2 ) be the probability to find a quark -antiquark pair inside a virtual photon of 
virtuality Q 2 and &(xg,r t ) be the total cross section for the dipole-nucleon interaction. 
The total cross section for the interaction of the virtual photon with the nucleus is given 
by iiii 



cT t ot(l*p)(x B ,Q 2 ) = J dzj d 2 r t P^(z,r t ;Q 2 )a(x B ,r t ) (2) 

P"\z,r t ;Q 2 ) = ^£^{l^r + |^| 2 } (3) 

/ 

\V T (z,r t -Q 2 )\ 2 = (z 2 + (l-z) 2 )a 2 K 2 1 (ar t )+m 2 f K 2 (ar t ) (4) 

\^ L (z,r t ;Q 2 )\ 2 = AQ 2 z 2 (l-z) 2 K 2 (ar t ) (5) 

where a 2 = z(l — z)Q 2 + m 2 . 

2. The optical theorem in the impact parameter b t representation reads 

a(x B ,r t )=2j d 2 b t lma e d l ipole (x B ,r t ,b t ) (6) 

where a el (x B ,b t ) is a partial elastic amplitude at fixed b t with a given Bjorken scaling 
variable x B . 

3. In all approaches to high energy scattering processes, it is assumed that the dipole - 
nucleus amplitude is purely imaginary at high energy, s-channel unitarity then leads to 
the following equation: 

/ n(x B ,r t ;b t ) \ 

^dipoie^B, t t - bt) = if 1 - e a (7) 



In the Glauber-Mueller model the opacity Q(x B ,r t ;bt) represents the exchange of one 
hard Pomeron. It is defined as 

7[ 2 QLg 4 

Vt(x B ,r t ;b t ) = K(x B ,r t ) S(b t ) nR 2 A = — — r 2 x B G(x B , ^) S{b t ) , (8) 

3 rf 
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so that for small Q(xB,r t ; b t ) Eq. (|6]) yields the DGLAP evolution equation (in the dou- 
ble ln(l/xs) ln(Q 2 )) approximation. We assume that this relation is valid in the whole 
kinematical region. 

4. We assume, for simplicity, a Gaussian form for the profile function 



1 

ttR 



i 



SiPt) = w a (9) 



where Ra is a target radius. 

5. Diffraction of the photon on the target can be thought of as an elastic scattering of each 
dipole off the nucleus So instead of lma dipole (xB, r t , b t ) we substitute 

\a d \ pole {x B ,r u b t )\ 2 = (l -e —J (10) 

to get (x D (7*A). 



It was shown in Ref. [|T6j, [I7|, |19[ that the eikonal approximation leads to a reasonable 
theoretical as well as experimental approximation in the region of not too small xb- We will 
use this model to calculate contributions of different twists to the total and diffractive cross 
sections. This needs additional justification which are based on the following arguments. 

1. When the shadowing corrections are small ( at large values of Q 2 ) the eikonal approx- 
imation gives the correct matching to the DGLAP evolution; 

2. The eikonal approximation gives a natural explanation and an estimate of the saturation 
scale Q fl (sb, A) [1,0,0; 



3. It is a good approximation |T^| to the non-linear evolution equation for the HERA kine- 
matical region at Q 2 > 5 GeV 2 where it describes available experimental data well. 

4. It provides the correct anomalous dimensions for all higher twist contributions in the 
limit N c > 1 and a s (Q 2 ) ln(l/x B ) > 1| 



5. Our approach includes the impact parameter behaviour of the dipole-target scattering 
amplitude and, therefore, it can be used for different target , including nuclei; 



6. The eikonal approach is very similar to the saturation model of Ref. JTD| which describes 
the experimental data at HERA very well. The main difference between these models is 
that eikonal satisfies the DGLAP evolution equation in the limit of large Q 2 , while the 
Golec-Biernat and Wusthoff model does not do so (for the complete discussion of this 
issue see Ref. 



This simple eikonal model is a reasonable starting point for understanding how the shadow- 
ing corrections are turned on. Moreover, as was shown in Ref. || using the Glauber approach 
for a nuclear target can be justified, there is a possibility that it will also be successful for the 
case of a nucleon (see Ref. [j2"T|). 
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2.2 Analytical expansion for mj = 0. 

Collecting formulae of the previous subsection we arrive at the following expressions (for three 
flavours) : 

atot(l*A)(x B ,Q 2 ) = 2^J dzj d 2 r t J d 2 b t {\^ T \ 2 +\^ L \ 2 } (l-e- 3 ^ 1 ) (11) 

a^A){x B ^Q 2 ) = ^JdzJ d\J d%{\* T \ 2 +\* L \ 2 } (l-e-^^) 2 (12) 

We introduce new dimensionless integration variables which simplify the calculation accord- 
ing to 

t 2 = r 2 Q 2 ] £ = e ^; y 2 = z(l - z) . (13) 
In the limit m 2 = Eq. (|TTD reads: 

* L (j*A) = Z2R\a em fJ 2 dyj™ dt 2 £ ^ -^L= K 2 (yt) (l - e~%) (14) 

* T {>f A) = m\a em fj 2 dy f dt 2 ^ | ^IZfJ, K i(yt) (l - e-T) (15) 
Vtot{l*A) = a L { 1 *A) + a T { 1 *A) (16) 

and a similar expression for Eq. (|I~2j). 

Most of the dipoles in the virtual photon wave function are of size r 2 = 4/Q 2 . We assume 
that all dipoles are of this size. This assumption will not affect our calculations of the twist 
contributions to the cross sections, as it only gives logarithmic Q 2 corrections at small Q 2 and 
is valid in the leading In Q 2 approximation. Thus, Eq. (§) reads 



n(x B ,t;£;Q 2 ) _t 2 £ 



where we have defined^ 



2 Q 2 

_ Aira s x B G(x B , Q 



C(x B ,Q 2 ) (17) 



C(x B ,Q)- 

6R A 

In this massless approximation integrals of Eq. ([n]) and Eq. flT2|) can be expressed in terms of 
special functions. 

First, integration over £ gives: 

_£ca (t 2 C\ „ (t 2 C\ ~ {-l) k+1 (t 2 C\ k , s 

where £i is an exponential integral and 7 is the Euler constant. We now consider the longitu- 
dinal part of the cross section. Integrating over t and summing over k one obtains 

dyy 5 ~ (-l) fc+1 / C\ fc 4T(l + k)\y 2 )- l - k 



vl{i*A) = y/^a em R 2 A f Vy ]T 

JO \/l — V 2 



16C a /■! rfyy 5 4C 



2a em ^R 2 A^Gl 2 _ ( Q-r 



em A Q 2 34 Uc 



0,1; I 



(20) 



' ' 2 

1,1,2;-^ 

1 In our model C(xb, Q 2 ) is a saturation scale. In general, the saturation scale has to follow from the solution 
of the complete non-linear evolution equation as we have already mentioned. 
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(a) (b) 

Figure 1: Integration contours for the Meijer function in Eq. (pOD and Eq. (p^). 



where we have used the definition of the generalized hypergeometric function 3^ and expressed 
the integral over y as a Meijer function G\\ (see Ref. |[22|| ). This Meijer function is defined by 
the following integration 



^32 
^34 



'91 

AC 



a li °-2] «3 



1 

2m Jl 



ds 



'g 2 V r(i- fll + s )r(i - a 2 + s )r(&i - s)r(b 2 - s)r(b 3 - s) 

AC 



T(a 3 - s)T(l -b 4 + s) 



(21) 

where the integration contour L runs from +00 to +00 and encloses all poles of functions 
r(6j— s), i — 1, 2, 3 in the negative direction, but does not encloses poles of functions r(l— cij+s), 
j = 1,2. The contour for the Meijer function in Eq. ( pO"P is shown in Fig.|l](a). 

In the case of the transverse cross section we repeat the by now familiar procedure to obtain 
the following result: 



C 

Aa em nR 2 A — 



n z2 ( Q 2 
Gs4 [AC 



M.i;i 



- -G 



3 2 
34 



AC 



0,1; § 

1,1,1;- 



+ 



AC f 1 
Qt\A 



G 



3 2 
34 



AC 



i.i; I 



2 2 2-1 

^, 2 



- -G 



3 2 
3 4 



AC 



1,2; J ' 
2,2,2; I. 



(22) 



The contour for the first two functions is shown in Fig.|l|(a) and for the last two in Fig.fjjfb). 

The simplest way to calculate cr^A)^ and a{^*A)^ is to note that using Eq. (0) and 
Eq. (0) 



d% 1 



_S2\2 

e 2 



e~) - / d\ (1-e 



^(2«) 



Using Eq. (|20|) and Eq. (|22|) one gets 



vtfAft^C/Q 4 ) = a(^A) TjL (C/Q') - -atf A) T , L {2C j Q>) 



(23) 



Formulae Eq. (|20[) , Eq. (|22|) and Eq. ( P3|) make it possible to estimate the contribution of 
any twist to the particular cross section. To this end we should close the contours of Fig.[TJ to 
the left over a particular number of poles. Each of these poles represents the contribution of a 
particular twist. Below we show some results for twist-2 and twist-4 contributions. 2 



a 



T=2 



7,2 C 8 
2a em nR A -- 



(24) 



2 Closing the contours to the right provides a small Q 2 series which, however is of little interest, since we 



assumed that 



0. 
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, , n 2 C 14 + 12( 7 + ln(QV4C)) 
a T = 4a em irR A — — , (25) 

°l t=2 = , (26) 



a^= 2 = Aa em 7rR 2 A —\n2 , (27) 



.I- = (|) 32(17 + 30 7 + 30 lnWV 4 C )) (2g) 

2.3 Ratios of cross sections 

As we now have analytical formulae for different twists we can estimate in which kinematical 
region twist-2, say, equals twist-4, and therefore find where the twist expansion breaks down. 
We will do this in the following section. The estimation we make is not accurate for two rea- 
sons. First, the equality of two first twists is a sufficient but not necessary condition for the 
breaking down of the twist expansion. This argument will be illustrated in the next section. 
Second, all models describing the shadowing corrections are written in the leading logarithmic 
approximation. In order to estimate the reliability of these models one has to perform calcu- 
lation of the corrections to the leading order. This has not been done yet. In Ref. |TB|, [H| we 



overcame this problem by noting that in the eikonal approximation corrections to the leading 
log factorize out of the formulae for the cross sections. This makes it possible to describe the 
experimental data quite well by introducing some universal factor. However, we would like to 
avoid introducing factors of unclear theoretical origin. 

One can overcome these difficulties by studying the ratios of the cross sections instead of 
the cross sections themselves. In this paper we will consider two ratios <Jl/&t and a® /a®- 
Both of these ratios have a remarkable property as a function of photon virtuality Q 2 . At small 
Q 2 the ratios vanish since the real photon is transverse. At large Q 2 they vanish as well, as 
implied by pQCD . This leads us to suggest that both ratios have maximum. We would expect 
that this maximum occurs at the same value of Q 2 max indicating that there exists a universal 
scale Q 2 oc Q 2 max which we call the "saturation scale". If in the kinematical region that we 
are interested in we get different maxima for different ratios, then we can conclude that the 
Shadowing Corrections are small in this specific kinematical region. The maximum is strongly 
influenced by the confinement forces at the scale m 2 . 

In practice it is difficult to separate the contributions of unitarity and confinement forces. 
The manifestation of this is that and erf are logarithmically divergent in the massless limit 
as Q 2 — > 0. We, however, expect that at large enough energies (small xb) the unitarity forces 
become important at much larger Q 2 then the confinement ones. This should appear as a 
maximum at virtualities much larger than the quark mass squared, which has been introduced 
to model confinement forces. 



3 Numerical estimations 

We now present the results of numerical calculations. We performed the calculation using 
GRV'94 parameterization for the gluon structure function xbG(xb, Q 2 )- The reason for choos- 
ing this particular parameterization is that we hope that non-linear corrections to the DGLAP 
evolution equation which were, perhaps, seen at HERA are not obscured by it. It is well known 
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that all existent experimental data can be fitted more or less well with a linear parameteri- 
zation. It seems senseless to look for shadowing corrections using such parameterization. An 
additional reason for using GRV'94 is that it enables us to perform calculations at values of Q 2 
as low as0 0.4 GeV 2 . 

However, GRV'94 has a serious drawback. There exists a domain in Q 2 where it gives 
the anomalous dimension of xbG(xb,Q 2 ) larger than 1. This manifests itself as a spurious 
minimum of 1/ k(xb,4/Q 2 ) (see Eq. ([8])) as a function of Q 2 at some Q 2 (xb)- This minimum 
threatens to spoil the calculation of extremum of ratios discussed in the previous section. To 
get rid of it we use the following gluon structure function 

x B G(x B ,Q 2 ) = x B G GRV (x B ,Q 2 )9(Q 2 -Q 2 (x B )) 

Q ~ x B G GRV (x B , Q\x B )) 6(Q 2 (x b ) - Q 2 ) (30) 



Q 2 (X B 



where 9 is a step function. Eq. ( p0|) is a rough way of taking into account the shadowing 
corrections in the gluon channel which have not been calculated systematically in our approach, 
namely, Eq. (|30"D leads to a DGLAP gluon structure function for Q 2 larger than the saturation 

scale Q 2 {xb) while it gives a saturated gluon structure function for Q 2 smaller than this scale. 

The results are conveniently represented in terms of various structure functions instead of 
corresponding cross sections. They are defined by the equation 

F 2 (x B ,Q 2 ) = -^a( 1 *A) 

Note, that all the equations of our model are written in the leading Iu(1/xb) approximation^. 
This means that they are valid in the kinematical region where the condition as{Q 2 ) \b.(1/xb) > 
1 holds. Hence, we are working at the boundary of this kinematical region. So the question 
concerning corrections to the eikonal approach arises. These corrections are large but, as we 
have discussed, they cancel in the ratios since they mostly affect the normalization of the 
calculated structure functions. We can take these corrections into account by modeling the 
anomalous dimension as described in Ref. [19]. This effectively leads to the factor 0.5 in 



front of Eq. ( p0|) and Eq. (|22|). To compare our calculation of the structure functions with 
experimental data one should divide all the values by factor 2. 

To estimate the contributions of the different twists to various structure functions we per- 
form calculations using Eq. @jj),Eq. fl2|) and Eq. fl2|). Note, that Eq. (H) and Eq. fl2|) are 



valid only in the massless limit m 2 f = 0. As was explained above, this is a good assumption as 



long as we are not concerned with small Q 2 behaviour. All our calculations are made for the 
eRHIC kinematical region. The results of the calculations for Zn (A=30), Sn (A=119) and U 
(A=238) are shown in Figs. |2|,0 and [|. The following important effects are clearly seen when 
we compare the contributions of first two non- vanishing twists: 

1. In all structure functions except Fl there is a scale Q 2 specific for the particular struc- 
ture function where these two twists are equal. Q 2 is largest for F®; 

2. Q 2 grows as xb decreases and A increases; 

3. twist-4 in Fl cancels numerically with twist-4 inF] F 2 in the region of Q 2 where twist-4 
is significant (low Q 2 ). This is a reason for measuring the polarized structure functions 
instead of the total one. 



3 See full discussion in Ref. fllq| . 

4 And leading In Q 2 as well, see explanations to the Eq. ( |l7| ) . 

5 This effect was also noted in Ref. || where another, though similar, modeljlij was used for SC. 
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In Figs. | and E we plotted the ratios F L /F T (x B , Q 2 , A) and F£ / (x B , Q 2 , A) corre- 
spondingly versus Q 2 . The calculations were performed with formulae Eq. ( |TTD and Eq. ([12|) 
not neglecting the mass of the quark rrif. As is expected on theoretical grounds]?! the ratio 
has a maximum which increases as x B decreases and A increases. It occurs at Qm ax 3> m 2 
which implies that we are sufficiently distant from the confinement region. However, it differs 
numerically for inclusive and diffractive ratios. This in turn implies that the SC will be still 
small at eRHIC. 

In Figs. [7j and || we show the behavior of Qln^ix-B-, A) as a function of x B and A. Note in 
Fig. | that Q 2 max ~ A - 25 . If the SC were absent we would expect that Q 2 max not to depend on 
A, while in saturation regime Q 2 max ~ A 2 / 3 . This confirms once more our conclusion that the 
kinematical region of eRHIC is intermediate between the linear and the saturation regimes. 



4 Conclusions 

The goal of our work is to understand to what extent the high parton density regime of QCD 
could manifest itself in DIS experiments at eRHIC where the density of partons could be large 
enough to make shadowing corrections significant despite the fact that energies are much smaller 
than that achieved at HERA for DIS on nucleons. 



The eikonal approach, developed in Refs. [12, 19, 15| was used for the first estimates of 



the manifestation of the saturation scale in deep inelastic scattering with nuclei. It has been 
argued that it provides a reasonable first iteration to the exact evolution equation. We realize 
that it is not clear whether the equations written in the leading ln(l/x B ) approximation are 
valid in the kinematical region of our interest. However, we believe that our calculations give 
a reasonable estimate for the magnitude of the effect and of the kinematic region where this 
effect could be observed experimentally. 

We observed two manifestations of the saturation scale: 

1. The position of maxima in ratios Fl/Ft and F® / F® (Qmox) ls hitimely related to the 
saturation scale Q 2 (x B ). Q 2 max shows x B and A dependence typical for the saturation 
scale (see Figs. 7 and 8); 

2. The results of our calculations show that there exists a scale at which the twist expansion 
breaks down since all twists become of the same order. We can use the DGLAP evolution 
equations only for Q 2 larger than this scale. The energy (x B ) and A dependence of this 
scale gives us hope that the experimental data for deep inelastic scattering with nuclei 
will help to separate leading and higher twists contributions. In the case of nucleon deep 
inelastic scattering such a separation appears to be a rather difficult task which has not 
been performed until now. 

With this paper we started a systematic analysis of possible manifestation of gluon satu- 
ration phenomena in deep inelastic scattering on nuclei. We hope that our estimates will be 
useful for planning of future experiments especially at eRHIC. 
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Fi gure 2: Different twists contributions to the various structure functions for DIS on nucleus A=30: leading 
twist ( at high Q 2 ) - dashed line, next-to-leading - dotted one, exact structure function - solid curve. 



11 



Fl(Q 2 ) 

x B = 1(T 2 x B = 5- 10~ 3 x B = 1(T 3 




2 4 6 8 10 2 4 6 8 10 2 4 6 8 10 12 14 



Fi gure 3: Different twists contributions to the various structure functions for DIS on nucleus A— 119: leading 
twist ( at high Q 2 ) - dashed line, next-to-leading - dotted one, exact structure function - solid curve. 
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Fi gure 4: Different twists contributions to the various structure functions for DIS on nucleus A=238 (U): 
leading twist ( at high Q 2 ) - dashed line, next-to-leading - dotted one, exact structure function - solid curve. 
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Figure 5: Ratio F L /F T versus Q 2 for different x B and A. 
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Figure 8: Scaling of Q 2 max with A: In of maxima of ratios of Fl/Ft and F® /F® versus In A for different xb- 
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